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Abstract 

A comparative study of fracture in Al is carried out by using quantum me- 
chanical and empirical atomistic description of atomic interaction at crack 
tip. The former is accomplished with the density functional theory (DFT) 
based Quasicontinuum method (QCDFT) and the latter with the original 
Quasicontinuum method (EAM-QC). Aside from quantitative differences, the 
two descriptions also yield qualitatively distinctive fracture behavior. While 
EAM-QC predicts a straight crack front and a micro-twinning at the crack 
tip, QCDFT finds a more rounded crack profile and the absence of twinning. 
Although many dislocations are emitted from the crack tip in EAM-QC, they 
all glide on a single slip plane. In contrast, only two dislocations are nucle- 
ated under the maximum load applied in QCDFT, and they glide on two 
adjacent slip planes. The electron charge density develops "sharp corners" 
at the crack tip in EAM-QC, while it is smoother in QCDFT. The physics 
underlying these differences is discussed. 

Key words: Plastic Deformation, Dislocations, First-Principles Electron 
Structure Theory, Atomistic Simulation, fracture mechanics 
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1. Introduction 

Understanding fracture behavior in materials is a challenging undertak- 
ing. Despite nearly a century of study, several important issues remain un- 
solved. For example, there is little fundamental understanding of brittle to 
ductile transition as a function of temperature in many materials; there is 
still no definitive explanation of how fracture stress is transmitted through 
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plastic zones at crack tips; and there is no complete understanding of the 
disagreement between theory and experiment regarding the limiting speed of 
crack propagation. These difficulties to a great extent stem from the fact that 
fracture phenomena are governed by processes occurring over a wide range of 
length and time scales; thes e processes are all connected and all co ntribute 
to the total fracture energy ( Van der Giessen and Needlemanl . 2002 ). 

As emphasized by Van de Giessen and Needleman, although the atomistic 
interaction at a crack tip may only account for a small fraction of the total 
fract ure energy, it can be a controlling factor (IVan der Giessen and Needlemanl . 
20021 ) - after all, all fractures take place by breaking atomic bonds. Critical 
atomistic information, such as surface energy, stacking fault energy, dislo- 
cation nucleation/propagation energy and twinning formation energy, etc, 
has been known to play central roles in fracture. In fact, some of these 
quantities are at the heart of fractu r e mec hanics, including t he Gr i ffith's cri- 
terion for brittle fracture ( Griffith . 1920 ). Rice's criterion (jRicel . 1992 ) for 
crack tip blunting and more recently a criterion for twinning at crack tip 
(ITadmor and Hail . 120031 ). to name a few. 

Because of the inherent multiscale nature of fracture - the process at 
each scale depends strongly on what happens at the other scales, the mod- 
eling and simula t ion o f fracture calls for concurrent multiscale approaches 
( Lu and Kaxirasl . 120051 ). One of the first concurrent multiscale modeling of 
fracture was based o n Macroscopic A t omist ic ab initio Dynamics (MAAD) 



method for Silicon (IBroughton et all Il999l ). MAAD couples a quantum 
mechanical description of atoms at crack tip, to an empirical (or classical) 
atomistic description of atoms at a short distance away from the crack tip, 
and to the continuum finite-element description of the rest of the system. 
Since MAAD, several other concurrent multiscale methods have been de- 
veloped, all involving s ome l evel of quantum mechanical modeling at th e 



crack tip (jCsanyi et all 12004 : iBernstein and Hessl . 120031 : lOgata et all l200ll ) . 



All these methods were developed/applied for Si owning to the following 
technical reasons: (1) large-scale electronic structure methods such as linear- 
scaling algorithms are only applicable to covalently-bonded semiconductors 
like Si; general approaches for metals still remain elusive; (2) satisfactory 
QM/M M coupling schemes f o r metals were l e ss well developed u ntil re- 
cently ( Bernstein etall . 120091 : iLiu et all 120071 : Izhang and Lul . l2007h . On 
the other hand, concurrent multisca le approache s that do not involves quan- 
tum mechanic s are readily available flMiller et all ll998l : lKohlhoff et allfl99~ll 



Buehler et all I2006T ). Among them, Quasicontinuum (QC) method is par- 
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ticularly promising and it has bee n widely applied to many materials prob- 
lems, including fracture in metals (ITadmor and Millerl . 120051 ) . QC strives to 
achieve a "seamless" coupling between atomistic and continuum descriptions 
and allows quantum mechanical interactions incorporated in a systematical 
manner. For example, although the original QC was based on classical atomic 
interactions, significant progress has been made rece ntly to incorporate quan- 
tum mechanical i nteractions in th e local QC region flHaves et all 120 06). non- 
local QC region (iLu et al. . 20061 ) and entire QC system ( Peng et al. . 2008 ). 
The coarse-graining strategy of QC ha s also been explored to perform large- 



scale electronic structure calculations (IGavini et all 120071 ) . 



Despite the impressive advance in multiscale methodology development, 
a crucial question remains unanswered. Although it is clear that an atomistic 
description at a crack tip is indispensable for many purposes, it is not well 
established whether a quantum mechanical description at the crack tip is 
truly necessary. This is a poignant point given the continuing improvement 
of empirical potentials and the still heavy costs for quantum simulations. It 
is to address this question that motivates the present study. As a first look 
at the problem, we focus on crack tip plasticity in Al and compare results 
received from a quantum mechanical description vs. an empirical atomistic 
description at the crack tip, both in the framework of QC. Since the atom- 
istic resolution is only necessary near the crack tip while the linear elastic 
fracture mechanics boundary conditions can be applied in the far field, QC is 
well poised for such fracture simulations. In addition, QC simulations involve 
quasi-static energy minimization, thus the unrealistically high strain rates 
common to molecular dynamics simulations are avoided. Unfortunately, as a 
result thermally activated processes are precluded in QC. Al is chosen in this 
study because it is relatively inexpensive for density functional theory (DFT) 
calculations and an excellent embedded-a t om-m ethod (EAM) empirical po- 
tential exists for Al (jErcolessi and Adamsl . 1 1994 ) . The goal of this work is to 
examine how and why the results received in the empirical simulations differ 
from those in quantum mechanical DFT simulations at the crack tip. To this 
end, we employ so-called QCDFT method in which the nonlocal atoms at 
the crack tip are treated with DFT. The QCDFT results are compared to 
those obtained from the original QC method in which all nonlocal atoms are 
treated with EAM empirical potential. 

The structure of the paper is as follows. The methodology is introduced 
in Section 2 for both QC and QCDFT methods. A semi-infinite crack un- 
der mode I loading is set up in Section 3.1. The computational parameters 
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are described in Section 3.2 and the loading procedure for the crack is sum- 
marized in Section 3.3. The simulation results and analysis are presented 
in Section 4 and discussions are given in Section 5. Finally we conclude in 
Section 6. 



2. Methodology 



The QC me thod (IShenov et all Il999l : iTadmor et all 1 19991 ) is a multi- 
scale approach ( Lu and Kaxirasl . 20051 ) that combines atomistic models with 



continuum theories, and thus offers an advantage over conventional atom- 
istic simulations in terms of computational efficiency. The idea underlying 
the QC method is that atomistic processes of interest often occur in very 
small spatial domains (e.g., crack tip) while the vast majority of atoms in 
the material behave according to well-established continuum theories. To 
exploit this fact, the QC method retains atomic resolution only where neces- 
sary and coarsens to a continuum finite element description elsewhere. This 
is achieved by replacing the full set of TV atoms with a small subset of N r 
"representative atoms" or repatoms (N r <C N) that approximate the total 
energy through appropriate weighting. The energies of individual repatoms 
are computed in two different ways depending on the deformation in their 
immediate vicinity. Atoms experiencing large deformation gradients on an 
atomic-scale are computed in the same way as in a standard fully-atomistic 
method. In QC these atoms are called nonlocal atoms. In contrast, the en- 
ergies of atoms experiencing a smooth deformation field on the atomic scale 
are computed based on the deformation gradient in their vicinity as befitting 
a continuum model. These atoms are called local atoms. The total energy 
E tot (which for a classical system can be written as E tot = ^f =1 E i) with Ei 
the energy of atom i) is approximated as 



N nl N lc 



= $>i({q}) + y £n j E^({F}). C 1 ) 

i=l j=l 

The total energy has been divided into two parts: an atomistic region of A^ nl 
nonlocal atoms and a continuum region of N loc local atoms (A^+A^ 100 = A^ r ). 

The original formulation of QC was limited to classical potentials for 
describing interactions between atoms. However, since many materials prop- 
erties depend crucially on the behavior of electrons, such as bond break- 
ing/forming at crack tips or defect cores, chemical reactions with impurities, 
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surface reactions and reconstructions, electron excitation and magnetism, etc, 
it is desirable to incorporate appropriate quantum mechanical descriptions 
into the QC formalism. QCDFT is one strategy to fill this role. In specific, 
QCDFT combines the coarse graining idea of QC and the coupling strategy 
of the quantum mechanics/molecular mechanics (QM/MM) method. This 
method can capture the electronic structure at the crack tip within the ac- 
curacy of DFT and at the same time reach the l ength-scale that is relevant 
to experiments (ILu et all 120061 : IPeng et all 120081 ). 

The original QC formulation assumes that the total energy can be writ- 
ten as a sum over individual atom energies. This condition is not satisfied 
by quantum mechanical models. To address this limitation, in the present 
QCDFT approach t he nonlocal regi o n is treated by a n EAM-based QM/MM 
coupling approach ( Lu et al. . 120061 : Liu et al. . 2007 ): the Kokn-Sham den- 
sity functional theory (KS-DFT) is coupled to EAM with the interaction 
energy calculated also by EAM. The local region, on the other hand, is dealt 
with by EAM, which is the same energy formulation used in the MM part 
of the nonlocal region. This makes the passage from the atomistic to contin- 
uum seamless since the same underlying material description is used in both. 
This description enables the model to adapt automatically to changing cir- 
cumstances (e.g. the nucleation of new defects or the migration of existing 
defects). The adaptability is one of its main strengths of QC and QCDFT, 
which is missing in many other multiscale methods. 

More specifically, in the present QCDFT approach the material of in- 
terest is partitioned into three distinct types of domains: (1) a nonlocal 
quantum mechanical DFT region (region I); (2) a nonlocal classical region 
where classical EAM potentials are used (region II); and (3) a local region 
(region III) that employs the same EAM potentials as region II. The cou- 
pling between regions II and III is achieved via the QC formulation, while the 
coupling between regions I and II is accomplished by the QM/MM scheme 
(ICholv et all l2005l:lLiu et all l2007h . The total energy of the QCDFT system 



is then (ILu et all . 120061 ) 



QCDFT _ ™i 



^ nl [i + n] + Eri^^] oc ({ F }) 



tot 



= -E/dft [I] - £eam[I] + £eam[I + II] + £f=T n,£j oc ({F}), 



(2) 



where E nl [l + II] is the total energy of the nonlocal region (I and II combined 
with the assumption that region I is embedded within region II), £"dft[I] 
is the energy of region I in the absence of region II computed with DFT, 
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£eam[H] is the energy of region II in the absence of region I computed with 
EAM, and Eeam[I + H] ^ s the energy of the nonlocal region computed with 
EAM. 

Other types of combination with quantum mechanical and classical atom- 
istic methods may also be implemented in QCDFT. The great advantage of 
the present implementation is its simplicity; it demands nothing beyond what 
is required for a DFT calculation and an EAM-QC calculation. Another im- 
portant practical advantage of QCDFT method is that, if region I contains 
many different atomic species while region II contains only one atom type, 
there is no need to develop reliable EAM potentials that can describe each 
species and their interactions. This is because if the various species of atoms 
are well within region I, the energy contributions of these atoms are canceled 
out in the total energy calculation. This advantage renders the method par- 
ticularly useful in dealing with impurities, which is an exceedingly difficult 
task for empirical potential simulations. 

The equilibrium structure of the system is obtained by minimizing the 
total energy in Eq. [2] with respect to all degrees of freedom. Because the 
time required to evaluate £?dft[I] is considerably more than that required 
for computation of the other EAM terms in E^ DFT , an alternate relaxation 
scheme turns out to be useful. The total system can be relaxed by using 
conjugate gradient approach on the DFT atoms alone, while fully relaxing 
the EAM atoms in region II and the displacement field in region III at each 
step. An auxiliary energy function can be defined as 

«}]= im ^[{q}], (3) 

{q n },{q m } 

which allows for the following relaxation scheme: (i) minimize E^ DFT with 
respect to the atoms in regions II ({q 11 }) and the atoms in region III ({q 111 }), 
while holding the atoms in region I fixed; (ii) calculate i? t ^f DFT [{q}], and the 
forces on the region I atoms; (iii) perform one step of conjugate gradient 
minimization of E'\ (iv) repeat until the system is relaxed. In this manner, 
the number of DFT calculations performed is greatly reduced, albeit at the 
expense of more EAM and local QC calculations. A number of tests have 
shown that the total number of DFT energy calculations for the relaxation 
of an entire system is about the same as that required for DFT relaxation of 
region I alone. 
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3. Computational details 

3.1. Model setup 




Figure 1: (Color online) (a) The overview of the entire crack system with finite-element 
mesh; (b) A blown-up view of (a) showing the nonlocal region; (c) region I box and atomic 
structure at the crack tip; (d) schematic partition of the system into region I, II and III. 
The x, y and z axis is along [111], [110], and [112], respectively. All lengths are in A. 



A semi-infinite crack in a single Al crystal is studied by both QC and 
QCDFT for comparisons. The crack is made by removing two layers of 
atoms with x < and y = 0&1.41 A. The crack plane is (110) and in 
this orientation, (111) plane is the only active slip plane for dislocations 
emitted from the crack tip; all other {lll}-type planes lie obliquely to the 
crack plane and are thu s precluded by the imposed plane strain conditions 
( ITadmor and Hail . 120031 ) . This configurat i on wa s used previously in a MD 
study by Hoagland et al. (jHoagland et all Il990l ) and an EAM-QC study by 
Hai et al. ( Hai and Tadmor . 20031 ) although the initial crack opening and the 
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EAM potential used are different in these studies. The initial crack opening 
in the present work is determined based on two competing considerations: (1) 
it cannot be too narrow otherwise the crack will close and/or large number 
of loading steps is required to observe the onset of plasticity; (2) it cannot 
be too wide otherwise the DFT region is too large to render calculations 
feasible. Of course, the DFT region has to be large enough to capture the 
crucial plasticity events at the tip. The crack is subject to mode I loading 
along y direction shown schematically in Fig. [H(d). 

The dimension of the system is 6000 x 6000 x 4.887 A 3 along the x, y, z 
directions respectively. The system is periodic in z-direction, and has Dirich- 
let boundary conditions in the other two directions. The system contains 
over 11 million Al atoms - a size that is well beyond the reach of any full- 
blown quantum calculation. The schematic overview of the system is shown 
in panel (a) of Fig. [1] with finite element meshes. Panel (b) is a zoomed-in 
view of the nonlocal region and panel (c) displays the atomic detail of the 
crack tip. 



3.2. Computational parameters 

Two comparative calculations are carried out for the crack: EAM-QC 
vs. QCDFT for Al. EAM-QC is the original QC with EAM potential for 
atomic interactions. In this work, the EAM potential used is res caled from 
the original "force-matching" EAM (jErcolessi and Adamsl . Il994l ) potential 
so that it matches precisely the value of the lattice c o nstan t and bulk mod- 



ulus of Al from the DFT calculations (jCholy et all 120051 ) . Although the 



re-scaling changes very little to the original potential, it eliminates lattice 
parameter mismatch at the QM/MM interface, and thus reduce QM/MM 
coupling errors. 

DFT calculations are carried out with the Vienna Ab-initio Simulation 



Package (VASP) (IKresse et allll993l : lKresse and Hafnerlll 994 : iKr esse and Furthullerl . 
Il996al lbl) which is based on Kohn-Shem Density Functional Theory (KS- 
DFT) with the local density approximation and ultrasoft pseudopotentials. 
A plane- wave cutoff of 129 eV is used in the calculat ions and the fc-points 
are sampled according to the Monkhorst-Pack method (Monkhorst and Pack, 



19761 ) with a 1 x 1 x 9 mesh in the Brillouin zone. There are 134 DFT atoms 
in region I. 
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3.3. Loading "procedure 

The simulations are performed quasi-statically with displacement bound- 
ary conditions where the displacement is prescribed as a function of intended 
stress intensity factor (SIF) at each loading step. For EAM-QC, the loading 
procedure is outlined as following: 

(1) For a small initial stress intensity f actor Kj, the anisotropic Linear Elastic 
Fracture Mechanics (LEFM) solution (jSih and Liebowitzl . ll968 ) Ulefm(X, Kj) 
is obtained. Each atom is displaced according to u(X) = Ulefm(X, Kj), 
where u is the displacement field and X is the position of nodes in the 
model. 

(2) The displacement at the model boundaries is fixed except for the crack 
surfaces; the positions for all repatoms are obtained by energy minimization 
as explained before. 

(3) The finite-element mesh, the status of repatoms (local vs. nonlocal) and 
the neighbor list are updated. 

(4) The SIF is increased by a small amount AK\ as K\ — K\ + AK\ and 
repeat from step (1) until the intended SIF is achieved. 



9 2 5 

In this study the increment AKi = 0.001 eV/A ' is u sed. The loading; 



proc edure adopted here follows closely that of reference (IHai and Tadmorl . 



20031 ) 



Because DFT calculations are much more expensive than EAM, we use 
EAM-QC to load the crack until an incipient plasticity is about to take place, 
at which point QCDFT is switched on. In other words, a QCDFT relaxation 
starts from a configuration that is obtained by EAM-QC for K\. QCDFT 
then increases the SIF by AK\. This is a reasonable approximation because 
EAM is known to give accurate results for deformations prior the onset of 
crack tip plasticity or the appearance of lattice defects. As will be shown 
later, since the critical load for the onset of plasticity from QCDFT is smaller 
than that from EAM-QC, the present loading strategy does not run the risk 
of missing incipient plasticity of QCDFT. 

4. Results and Analysis 

4.1. EAM-QC calculation of Al 

To establish the validity of present EAM-QC method, we first perform 



EAM -QC calculations for the same crack studied in reference ( IHai and Tadmorl . 



20031 ) . We use the same EAM potential and the same initial crack opening as 
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Figure 2: EAM-QC results at SIF Kj = 0.144 eV/A 2 5 . (a)The out-of-plane displacement 
U z ; the displacement contours range from -0.5 A(darkest) to 0.5 A(lightest). (b)The 
zoomed-in view of the crack tip atomic structure and finite-element mesh. The open circle 
represents the atomic position. All distances are in A. 



that in (IHai and Tadmorl . 120031 ) . The results are presented in Fig. 14.11 where 
two edge dislocations are emitted from the crack tip and they glide at the 
same slip plane in a symmetrical manner. The critical SIF is 0.144 eV/ A 2 5 . 
All these results are identical to those found in (IHai and Tadmorl . 120031 ) and 
thus validate the present EAM-QC method. 

Next, we apply EAM-QC to the crack of interest. The crack is loaded 
quasi-statically and no crack tip plasticity is observed until the SIF reaches 

9 2.5 

K\ = 0.180 eV/A ' . By comparison, the critical SIF for pure brittle cleavage 

o 2 5 

computed from the Griffith criterion for this orientation is K\ c = 0.205 eV/A ' . 
In Fig. [3l we present the out-of-plane displacement U z as a function of ap- 

9 2 5 

plied K\ values. For K\ = 0.179 eV/A ' , although no dislocation is ob- 
served, significant deformation at the crack tip is clearly visible (panel a). 

o 2 5 

At K\ = 0.180 eV/A ' , the first dislocation is nucleated and subsequently 
moves away from the crack tip. The dissociated 1/2 [110] edge dislocation is 
stabilized at about 70 A below the crack on a (111) plane (panel b). The 
contour shading of Fig. [3] corresponds to the magnitude of U Z) whose non- 
zero values indicate that the edge dislocation is dissociated into partials. 
Interestingly, at the same time, a micro- twin is also nucleated from the crack 
tip (Fig. [5K) - The micro- twn is two layers in both length and width, the 
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Figure 3: The out-of-plane displacement U z obtained from the EAM-QC calculations at 
SIF Ki of (a) 0.179, (b) 0.180, (c) 0.184, (d) 0.186, (e) 0.191, and (f) 0.198 eV/A 2 " 5 
respectively. The displacement contours range from -0.5 A(darkest) to 0.5 A(lightest). 
The crack surfaces are represented by red curves. All distances are in A. 
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twinning direction is [714] and the twinning plane is (131) . As K\ value is 
increased, more dislocations are emitted on the (111) plane and at the same 
time, the micro-twin grows in length but not in width. More specifically, as 
Ki increases to 0.184, 0.186, 0.191, 0.198 eV/A 2 ' 5 , two, three, four and five 
dislocations are emitted from the crack tip and they glide on the same (111) 
plane, as shown in the panel of (c), (d), (e) and (f) of Fig. [3] respectively. 
Correspondingly, the micro-twin grows to three, four, five and six layers in 
length, respectively. The micro-twin structures of two, three, five and six lay- 
ers in length are shown in the panel of (a), (b), (c), (d) of Fig. [5] respectively. 
The width of the micro-twin is not increasing perhaps due to unfavorable 
stacking fault energy along the usual twin ning; plane. The reaso n that the 
twinning was not observed in Hai et al. ( IHai and Tadmorl . 120031 ) is proba- 
bly due to the the different initial crack openings rather than the different 
EAM potentials used. We have done additional calculations for the present 
crack opening (2 layers ) with the same EAM potential used in reference 
(IHai and Tadmorl . 120031 ) and found a similar twinning at the crack tip. The 
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Figure 4: The dislocation distribution function D(x) as a function of the inverse of the 
distance x. The dashed curve is a fit to the filled circles. 

emitted dislocations share the following characteristics: (1) They are all edge 
dislocations of 1/2[110] Burgers vector and dissociated into Shockly partials 
with a separa tion distance of 16 A. Th is result agrees with the previous study 
of Hai et ah ( IHai and Tadmorl . 120031 ) with the similar EAM potential. But 
this value is too large compared to the e xperimental splitting d i stanc e of 5.5 
A, measured by Mills and St adelmannf Mills and Stadelmannl . 1 19891 ). The 
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discrepancy may be attributed to still too low stacking fault energy of the 
EAM potential. (2) They are all on the same {111} slip plane whose position 
is x = a/2, where is a is (111) inter-plane distance. The active slip plane is 
slightly ahead of the crack front position at x = 0. The emitted dislocations 
move away from the crack tip and form a pile-up against the local/nonlocal 

o 2 5 

QC interface. For K\ = 0.198 eV/A ' , the dislocation density D(x) defined 
as the number of dislocations per unit distance al ong pileup line, is foun d 
to be the square root of the inverse of the distance Mirth and Lothelll970h . 
The distance x refers to the distance between the dislocation center and the 
local/nonlocal interface. The fitted curve (dashed line) in Fig. 0] qualita- 
tively agrees well with the elastic theory where the fitting function is given 
as D(x) = 0.13236(l/x) 1 / 2 . 
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Figure 5: Atomic structure at the crack tip from EAM-QC. Some of atomic planes are 
highlighted in pink to indicate the twin. Dashed lines represent the two-layer twin, (a) K\ 
= 0.180, (b) Kj = 0.184, (c) Kj = 0.191, (d) Kj = 0.198 eV/A 2 * 5 . All distances are in A. 
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4-2. QCDFT calculation of Al 

o 2.5 

The crack is first loaded by EAM-QC until K\ reaches 0.169 eV/A ' at 
which point QCDFT is started. This critical loading is determined by trial 
and error; we launch a number of QCDFT calculations at different K\ from 
previously relaxed EAM-QC structures and examine whether any crack tip 
plasticity is taking palce. The smallest K\ value that results in incipient 
plasticity is the critical loading. In comparison, the critical SIF for pure 

o 2 5 

brittle cleavage computed from the Griffith criterion is K\ c = 0.267 eV/A ' . 




-0.5 -0.3 -0.1 0.1 0.3 0.5 



Figure 6: The out-of-plane displacement U z obtained from QCDFT calculations at SIF 
Ki of (a) 0.169, (b) 0.170, (c) 0.174, (d) 0.175 (e) 0.178 eV/A 2 ' 5 respectively. The dis- 
placement contours range from -0.5 A(darkest) to 0.5 A(lightest). All distances are in 
A. 

The crack tip behavior of QCDFT is much different. At K\ = 0.169 

9 2 5 

eV/A ' , two dissociated edge dislocations are nucleated - one above the 
crack plane and one below it. In contrast to EAM-QC, the two dislocations 
glide at two adjacent (111) slip planes, as shown in Fig. 7(d). The positions 
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of the two slip planes are at x = —a/2 and x = — 3a/2 respectively; they 
are slightly behind the crack front position. The separation distance of the 
two Schockly partials is 16 A, the same as that in EAM-QC calculation. 
This result is consistent with the fact that the dislocation cores are outside 
the DFT box where atomic interaction is determined by the rescaled EAM 
potential. 

As the load is increased, the two emitted dislocations are driven further 
away from the crack plane, however, no more dislocation is nucleated within 
the maximum load that we have explored in this study. For the largest 
K\ = 0.178eV/A ' , there are only two emitted dislocations gliding at 380 A 
(up) and 258 A (down) away from the crack plane. Due to the computational 
cost of QCDFT, we did not pursue more calculations for even larger loadings. 
However, it is evident from the present study that the crack tip plasticity 
observed in QCDFT is qualitatively different from that in EAM-QC. The 
differences are more striking given the fact that the QCDFT results are 
continued relaxation of EAM-QC configurations. 

5. Discussion 



Table 1: Relevant quantities calculated by VASP and EAM for bulk Al; the corresponding 
experimental values are extrapolated to T=0 K. 
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5.1. Dislocation nucleation at the crack tip 

EAM-QC shows that the crack blunts by emitting dislocations gliding 
on a single slip plane, and the number of emitted dislocations increases as 
the SIF is increased. However, QCDFT predicts that two dislocations are 
nucleated from the crack tip and they glide at two adjacent slip planes. 
Moreover, the number of emitted dislocations remains the same (two) for all 
SIFs. These distinctions are highlighted in Fig. [7| where schematic diagrams, 
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Figure 7: The schematic diagram, displacement contours and atomic structure with finite- 
elment mesh showing the dislocation pattern at the crack tip. (a), (b) and (c) are the 
schematic diagram, displacement contour plot and atomic structure respectively from 
EAM-QC calculation for Kj = 0.198 eV/A 2 ' 5 . (d), (e) and (f) are the same from QCDFT 
calculation for Kj = 0.169 eV/ A 2 ' 5 . 



16 



displacement contours and atomic structures at the crack tip for the two 
calculations are presented. In the schematic diagrams (a and d), the crack 
tip is represented by solid black lines/curves, and the dotted black lines 
denote the three relevant (111) slip planes at x — —3/ 2a, —a/2 and a/2. The 
stacking fault between the Shockley partials is indicated by a short red line 
segment. The displacement contour plots are the same as those shown in 
Fig. [3] and Fig. [6j they are reproduced here for convenience to the reader. 
The open circles in the atomic structure plots represent atomic positions. 
The sheared finite-element mesh in the atomic structure plots is the result of 
passing-by dislocations (the amount of shear corresponds to the net Burgers 
vector of the dislocations). The crack tip profile is approximated by blue line 
segments in the atomic structure plots. In Fig. 7(c) the crack tip profile is 
approximated by rectangular line segments AB+BD+DE to model the fact 
that the crack tip is blunted by emitting dislocations on a single slip plane. 
The active slip plane is represented BD whose position is at x = a/2. In Fig. 
7(f), however, the crack tip profile of QCDFT is approximated by zigzag line 
segments AC + CE because two adjacent slip planes are activated, and their 
positions are at x = —a/2 and x = — 3 /2a. 




Figure 8: A simple model for crack tip profile. The dotted line represents the relevant slip 
planes near the crack tip. The rectangular line segments AB+BD+DE and the zigzag 
segments AC + CE approximate EAM-QC and QCDFT crack tip profile respectively. 

To understand the origin of the differences, we resort to a simple model 
that captures the essential features of the crack configuration shown in Fig. El 



17 



One can consult Fig. [3 to understand the correspondence between the model 
and the actual crack tip configuration. The reason why the rectangular seg- 
ments AB+BD+DE are preferred in EAM-QC while the zigzag segments 
AC+CE are favored in QCDFT can be understood from the following sur- 
face energy analysis. Without loss of generality, we consider here the case 
where two dislocations are emitted for both EAM-QC and QCDFT. The 
length of BC equals to the magnitude of the Burgers vector because a full 
dislocation has been emitted upper- ward from A. As a result, AC repre- 
sents a {210} surface. There are two layers of atoms removed in the initial 
crack openning, which is one Burgers vector wide. Therefore in addition 
to a full dislocation emitted downward from E, the length of CD equals to 
twice of the Burgers vector magnitude. Hence CE is also a {210} surface. 
Therefore ~BD, AC (and CE) and AB (and KD) represents (111), {210} 
and (110) surface, respectively. The total surface energy associated with the 
zigzag and rectangular segments can be written as (AC + CE)j 2 io^o and 
(t4£>7ho + BD '7m + DEju )zo, respectively. z is the repeat distance along 
z axis. The various surface energies have been calculated by both DFT and 
rescaled EAM as summarized in Table 1. For this particular crack configu- 
ration, the surface energy of the rectangular segments is 0.05 eV lower than 
that of the zigzag segments based on EAM energetics. On the other hand, the 
surface energy of the zigzag segments is 1.23 eV lower than that of the rect- 
angular segments according to the DFT energetics (see note0). Therefore, 
the zigzag segments are preferred in QCDFT while the rectangular segments 
are favored in EAM-QC. The same conclusion holds for other loadings as 
well. 

5.2. Deformation twinning at the crack tip 

According to the Pe ierls criterion of deformation twinning at a crack tip 



( jTadmor and Hail . 120031 ). one can define twinnability which is the likelihood 



of a material to twin as opposed to slip at the crack tip. The d imensionless 



twinnability can be expressed as (jTadmor and Bernsteinl . 120041 ) 




T t = [1.136-0.151^]./^ (4) 

7us 



x a = y/l/3ao where ao is the lattice constant. AB = a, BC = ^3/2a, BC = y/Ea, 
~DE = 2a, AC = y/hj2a, ~CE vTOa, and z = y/§]2a. Taking a 2.3036 A, the total 
surface energy of the zigzag configuration is 3.67 eV for EAM and 3.93 eV for DFT. The 
total surface energy of the rectangular segments is 3.62 eV for EAM and 5.16 eV for DFT. 
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The coefficients 1.136 and 0.151 are universal constants for an fee lattice. 7 s f , 
7 US , and 7 ut are intrinsic stacking fault, unstable stacking fault and unsta- 
ble twinning energy respectively. A material will emit a dislocation before 
twinning if r t < 1 and will twin first if r t > 1. Our DFT and EAM calcula- 
tions find that r t is less than 1 as shown in Table I, which suggests that no 
true deformation twinning along (111) plane be formed, consistent with the 
QCDFT and EAM-QC simulation results. However, a two-layer-wide micro- 
twin along (131) plane is formed in EAM-QC while no such twin is present 
in QCDFT. The distinction between the QCDFT and EAM-QC results is 
due to the large discrepancy of 7 7f , o, wh i ch is t he energy barrier for a leading 
partial nucleation at a crack tip (jRicel 1 19921 ). The large DFT value of 7 US 
renders the nucleation of the leading partial difficult, which in turn prohibits 
the formation of deformation twinning. 




Figure 9: The electron charge density (A 3 ) at the crack tip for (a) EAM-QC at K\ = 
0.198 eV/A 2 ' 5 and (b) QCDFT at Ki = 0.169 eV/A 2 * 5 . The density contours range from 
(blue) to 0.26 (red). The blue sphere stands for atomic position and dashed line is an 
approximate profile of the charge density contours. 
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5.3. Electron density at the crack tip 

In Fig. 10, we present the electron charge density around the crack tip 
for the EAM-QC configuration of K Y = 0.198 eV/A 2 ' 5 (left) and the QCDFT 

o 2.5 

configuration of K\ = 0.169 eV/A ' (right). The charge density for EAM- 
QC is determined by a superposition of atomic densities (obtained by VASP) 
centered at each EAM atoms. The distortion of charge density due to the 
defects is clearly visible. In EAM-QC the atomic bonding is weakened along 
the twin plane while in QCDFT the atomic bonding is significantly disrupted 
at the center of the crack. More importantly, the charge density profile 
of QCDFT is smoother than that of EAM-QC as indicated by the dashed 
curves. The "sharp" corners of EAM-QC charge density lead to higher kinetic 
energy of electrons. Since EAM-QC does not involve quantum mechanics, the 
"sharp corner" is not penalized energetically and thus permissible. Of course, 
the electron charge density profile reflects the underlying atomic structure: 
"sharp corners" correspond to a straight crack front thanks to a single active 
slip plane; smooth corners correspond to a more rounded crack front. 

6. CONCLUSION 

In summary, we have carried out a comparative study of fracture in Al by 
using two distinctive atomic interactions: quantum mechanical density func- 
tional theory and empirical embedded atom method. The DFT description 
of the crack tip is achieved by QCDFT method while the empirical descrip- 
tion by EAM-QC method. In addition to quantitative differences, qualita- 
tively different fracture behavior is also observed between the two methods. 
EAM-QC predicts a more or less rectangular crack tip configuration while 
QCDFT yields a more rounded tip profile. The difference is due to the fact 
that the emitted dislocations glide on a single slip plane in EAM-QC while 
two adjacent slip planes are active in QCDFT. As the stress intensity factor 
is increased, more and more dislocations are emitted from the crack tip in 
EAM-QC while the number of dislocations remains the same up to the max- 
imum loading applied in QCDFT calculations. A micro-twin is observed at 
the crack tip in EAM-QC, but it is absent in QCDFT. The electron density 
profile at the crack tip is also different between EAM-QC and QCDFT. All 
these differences can be understood in terms of defect energetics, including 
surface energy and stacking fault energy. 

The different results received suggest that the atomic nature of a crack 
tip is important and an accurate description of the atomic interaction at the 
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crack tip is indispensable. Although empirical potentials can be developed 
by fitting to DFT results, it is unlikely they will reproduce all relevant DFT 
energetics. This is particularly so since one does not know a priori what 
are the relevant energetics for a given crack. If several chemical species are 
present in a crack tip, the task of fitting a satisfactory potential becomes 
even more daunting. Therefore the solution lies at an explicit quantum me- 
chanical description of the crack tip, most likely in a form of DFT-based 
multiscale modeling, such as QCDFT. The present paper concerns atomistic 
aspect of fracture which is important for many purposes. However, there are 
interesting fracture phenomena that do not depend on atomistic features and 
thus are beyond the scope of the present paper. Finally, we have not touched 
upon fracture dynamics. The questions - such as will finite temperature dy- 
namics amplify or diminish the differences that we observed and what are 
the best strategies to implement dynamics in a multiscale setting - remain 
unanswered. We hope that the present paper could spawn more research 
effort in answering these questions. 
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